function [patch] = CreateCircularGratingPatch(diameter, frequency, angle, phase, falloff, bbox, offset)
%  CreateCircularGratingPatch() - Create a 2D image of a circular grating
%
%    [patch] = CreateCircularGratingPatch2(diameter, angle, frequency, phase, falloff, bbox, offset)
%
%MANDATORY INPUTS:
%    diameter - grating diameter (in pixels)
%   frequency - grating frequency (in cycles per pixel, meaningful only when <.25)
%OPTIONAL INPUTS:
%       angle - grating angle (in degrees, clockwise, default: 0)
%       phase - grating phase (in cycles, default: 0)
%     falloff - 90/10% sigmoidal falloff distance (in pixels)
%        bbox - bounding box size ([x y] in pixels)
%      offset - offset from bounding box center ([x y] in pixels)
%OUPUT:
%       patch - grating patch

%Changelog
%  [2008/02/01] Valentin Wyart (valentin.wyart@chups.jussieu.fr)
%  [2011/03/15] K. N'DIAYE (kndiaye01@yahoo.fr)

if nargin < 7 || isempty(offset)
    offset = [0 0];
end
if nargin < 6
    bbox = [];
end
if nargin < 5 || isempty(falloff)
    falloff = 2;
end
if nargin < 4 || isempty(phase)
    phase = 0;
end
if nargin < 3 || isempty(angle)
    angle = 0;
end
if nargin < 2
    error('Not enough input arguments.');
end
switch numel(frequency)
    case 1
        stripes_lambda = log(1/0.01 - 1)/sin(2*pi*frequency);
    case 2
        stripes_lambda = log(1/0.01 - 1)/sin(2*pi*frequency(2));
        frequency = frequency(1);
    otherwise
        error('CreateCircularGratingPatch:NoFrequency','Spatial frequency must be provided')
end

diameter = floor(diameter/2)*2;
radius = floor(diameter/2);

if frequency > .25
    warning('CreateCircularGratingPatch:SaturatingFrequency','Frequency above 1/4 cycle per pixel are meaningless')
    %frequency = 0.25;
end

angle = mod(angle, 180);
% Compute angle from vertical, clockwisely
angle2 = angle+90;

phase = mod(phase, 1);

falloff = max(falloff, 1);

if isempty(bbox)
    bbox = [diameter diameter] + falloff + 2*abs(offset);
end
bbox = floor(bbox/2)*2;
cbox = floor(bbox/2) + 0.5;

[x y] = meshgrid((1:bbox(1)) - cbox(1) - offset(1), (1:bbox(2)) - cbox(2) - offset(2));

patch = zeros([bbox([2 1]) 2]);

% the grating
patch(:,:,1) = cos(2*pi*(frequency*(-sin(pi/180*angle2)*x + cos(pi/180*angle2)*y) + phase));
% rescale between 0 and 1
% patch(:,:,1) = patch(:,:,1)/2+.5;
% smooth through sigmoid

stripes_lambda = log(1/0.01 - 1)/sin(2*pi*frequency);
patch(:,:,1) = sigmoid(stripes_lambda*patch(:,:,1), [0 1]);

aperture_lambda = log(1/0.01 - 1)*2/falloff;
patch(:,:,2) = sigmoid(aperture_lambda*(sqrt(x.^2 + y.^2) - radius), [1 0]);

patch = round(patch*255);

end

function [y] = sigmoid(x, lims)
y = diff(lims)./(1 + exp(-x)) + lims(1);
end
